##################################################
# Replication Code
# Taeyong Park and Andrew Reeves
# "Local Unemployment and Voting for President: Uncovering Causal Mechanisms"
# Summary: Replicate Figure 4 in the paper.
##################################################

rm(list = ls())
library(foreign)
source("intervalEst.R") # This function produces a summary table of the estimated results 


# Effect of decreasing median household income by $1,000
# 2008 
dataInc08=read.table("output_inc08")
results08_inc=intervalEst(-dataInc08) # Negative sign to compute effects of decreasing income

# 2012 
dataInc12=read.table("output_inc12")
results12_inc=intervalEst(-dataInc12) # Negative sign to compute effects of decreasing income


# Effect of increasing gas prices by 10 cents
# 2008 
dataGas08=read.table("output_gas08")
results08_gas=intervalEst(dataGas08)
results08_gas[[1]]=results08_gas[[1]]/10  
results08_gas[[2]]=results08_gas[[2]]/10

# 2012 
dataGas12=read.table("output_gas12")
results12_gas=intervalEst(dataGas12)
results12_gas[[1]]=results12_gas[[1]]/10
results12_gas[[2]]=results12_gas[[2]]/10



#################################
# PLOT
#################################

ruler<-3:1
nRuler<-length(ruler)

par(mfrow=c(1,2), 
    mar=c(2, 0.5, 5, 0.5), 
    oma = c(2, 2, 0, .1), cex=1.1)
plot(x=results08_inc[[1]], y=ruler, 
     xlab="",
     ylab="",
     main="2008",
     xlim=c(-.005, .005), 
     ylim=c(0.75, 3.25),
     pch=16, cex=0.9,
     axes=F)
for(i in 1:nRuler){
  segments(x0=results08_inc[[2]][(i-1)*2+1],
           x1=results08_inc[[2]][(i-1)*2+2],
           y0=ruler[i],
           y1=ruler[i])
}
axis(3, at=c(-.005, 0, .005),
     labels=c("-0.5%", "0", "0.5"))
axis(1, at=c(-.005, 0, .005),
     labels=c("-0.5%", "0", "0.5"))
abline(v=0, lty=2)

plot(x=results12_inc[[1]],
     y=ruler,
     xlab="",
     ylab="",
     main="2012",
     xlim=c(-.005, .005), 
     ylim=c(0.75, 3.25),
     pch=16, cex=0.9,
     axes=F)
axis(3, at=c(-.005, 0, .005),
     labels=c("-0.5%", "0", "0.5"))
axis(1, at=c(-.005, 0, .005),
     labels=c("-0.5%", "0", "0.5"))
abline(v=0, lty=2)
for(i in 1:nRuler){
  segments(x0=results12_inc[[2]][(i-1)*2+1],
           x1=results12_inc[[2]][(i-1)*2+2],
           y0=ruler[i],
           y1=ruler[i])
}

mtext("AME", side=2, adj=0.68, outer=T, cex=1.2)
mtext("ADE", side=2, adj=0.42, outer=T, cex=1.2)
mtext("ATE", side=2, adj=0.17, outer=T, cex=1.2)
mtext("The likelihood of voting for the incumbent party", 
      side=1, padj=1, outer=TRUE, cex=1.2)
mtext("(A) Effect of Decreasing Median Household Income", 
      side=3, padj=2, outer=T, cex=1.5)



#################################

# Mediator: GAS #

#################################

ruler<-3:1
nRuler<-length(ruler)

par(mfrow=c(1,2), 
    mar=c(2, 0.5, 5, 0.5), 
    oma = c(2, 2, 0, .1), cex=1.1)
plot(x=results08_gas[[1]], y=ruler, 
     xlab="",
     ylab="",
     main="2008",
     xlim=c(-.01, .01), 
     ylim=c(0.75, 3.25),
     pch=16, cex=0.9,
     axes=F)
for(i in 1:nRuler){
  segments(x0=results08_gas[[2]][(i-1)*2+1],
           x1=results08_gas[[2]][(i-1)*2+2],
           y0=ruler[i],
           y1=ruler[i])
}
axis(3, at=c(-.01, 0, .01),
     labels=c("-1%", "0", "1"))
axis(1, at=c(-.01, 0, .01),
     labels=c("-1%", "0", "1"))
abline(v=0, lty=2)

plot(x=results12_gas[[1]],
     y=ruler,
     xlab="",
     ylab="",
     main="2012",
     xlim=c(-.01, .01), 
     ylim=c(0.75, 3.25),
     pch=16, cex=0.9,
     axes=F)
axis(3, at=c(-.01, 0, .01),
     labels=c("-1%", "0", "1"))
axis(1, at=c(-.01, 0, .01),
     labels=c("-1%", "0", "1"))
abline(v=0, lty=2)
for(i in 1:nRuler){
  segments(x0=results12_gas[[2]][(i-1)*2+1],
           x1=results12_gas[[2]][(i-1)*2+2],
           y0=ruler[i],
           y1=ruler[i])
}

mtext("AME", side=2, adj=0.68, outer=T, cex=1.2)
mtext("ADE", side=2, adj=0.42, outer=T, cex=1.2)
mtext("ATE", side=2, adj=0.17, outer=T, cex=1.2)
mtext("The likelihood of voting for the incumbent party", 
      side=1, padj=1, outer=TRUE, cex=1.2)
mtext("(B) Effect of Increasing Gas Prices", 
      side=3, padj=2, outer=T, cex=1.5)

